/*
 * get_potential.c
 *
 *  Created on: Apr 15, 2010
 *      Author: marchi
 */
#include "stdio.h"
#include "statutil.h"
#include "typedefs.h"
#include "math.h"
#include "vec.h"

#include "stdlib.h"
#include "get_potential.h"

int macro_filter(int nx, int ny, int nz, matrix box, real cut){
	int rank=DIM;
	int n0[DIM];

	n0[0]=nx;
	n0[1]=ny;
	n0[2]=nz;

	real dx=box[XX][XX]/(real) nx;
	real dy=box[YY][YY]/(real) ny;
	real dz=box[ZZ][ZZ]/(real) nz;

	flter_k=(fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*nx*ny*nz);
	flter_r=(fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*nx*ny*nz);

	fftwf_plan forward =fftwf_plan_dft(rank,n0,flter_r,flter_k,FFTW_FORWARD, FFTW_ESTIMATE);

	int mx= (int) rintf(0.5*cut/dx);
	int my= (int) rintf(0.5*cut/dy);
	int mz= (int) rintf(0.5*cut/dz);
	real fact=1.0/(real) ((2*mx+1)*(2*my+1)*(2*mz+1));

	for(int i=-mx;i<=mx;i++){
		int ia=(i>=0)?i:i+nx;
		for(int j=-mx;j<=mx;j++){
			int ja=(j>=0)?j:j+ny;
			for(int k=-mx;k<=mx;k++){
				int ka=(k>=0)?k:k+nz;
				flter_r[ka+nz*(ja+ny*ia)][0]=fact;
				flter_r[ka+nz*(ja+ny*ia)][1]=0.0;
			}
		}
	}
	fftwf_execute(forward);
	for(int i=0;i<nx*ny*nz;i++){
		flter_k[i][1]=-flter_k[i][1];
	}
	return 0;
}

int get_potential(real ***rho_in, matrix oc, int nx,int ny,int nz, real Volume, real ***phi_out, bool eFilter){

	int rank=DIM;
	int n0[DIM];
	int nfx,nfy,nfz;
	nfx=(nx % 2 == 0)? nx/2: nx/2+1;
	nfy=(ny % 2 == 0)? ny/2: ny/2+1;
	nfz=(nz % 2 == 0)? nz/2: nz/2+1;

	n0[0]=nx;
	n0[1]=ny;
	n0[2]=nz;

	fftwf_complex* rho;
	fftwf_complex* rhok;
	fftwf_complex* phi;

	rho =(fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*nx*ny*nz);
	rhok=(fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*nx*ny*nz);
	phi =(fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex)*nx*ny*nz);

	fftwf_plan forward =fftwf_plan_dft(rank,n0,rho,rhok,FFTW_FORWARD, FFTW_ESTIMATE);
	fftwf_plan backward=fftwf_plan_dft(rank,n0,rhok,phi,FFTW_BACKWARD, FFTW_ESTIMATE);


	real* r;
	r=&rho_in[0][0][0];
//	FILE *fp;
//	fp=fopen("pippo.dat","w+");
	int i,j,k;
	for(i=0;i<nx*ny*nz;i++){
		rho[i][0]=r[i];
		rho[i][1]=0.0;
		phi[i][0]=0.0;
		phi[i][1]=0.0;
	}
	real mw1,mw2,mw3,mw,fact;

	fftwf_execute(forward);

	int ia,ja,ka;
	for(i=0;i<nx;i++){
		ia=(i<nfx)?i : i-nx;
		for(j=0;j<ny;j++){
			ja=(j<nfy)?j : j-ny;
			for(k=0;k<nz;k++){
				ka=(k<nfz)?k : k-nz;
				int NN=k+nz*(j+ny*i);
				mw1=oc[XX][XX]*ia+oc[XX][YY]*ja+oc[XX][ZZ]*ka;
				mw2=oc[YY][XX]*ia+oc[YY][YY]*ja+oc[YY][ZZ]*ka;
				mw3=oc[ZZ][XX]*ia+oc[ZZ][YY]*ja+oc[ZZ][ZZ]*ka;
				mw=sqr(mw1)+sqr(mw2)+sqr(mw3);
				if(!j && !j && !k) fact=0.0;
				else fact=1.0/(mw*Volume*M_PI);
				rhok[NN][0]=fact*rhok[NN][0];
				rhok[NN][1]=fact*rhok[NN][1];

//				rhok[k+nz*(j+ny*i)][0]=fact*rhok[k+nz*(j+ny*i)][0];
//				rhok[k+nz*(j+ny*i)][1]=fact*rhok[k+nz*(j+ny*i)][1];

				if(eFilter) {
					rhok[NN][0]=rhok[NN][0]*flter_k[NN][0]-rhok[NN][1]*flter_k[NN][1];
					rhok[NN][1]=rhok[NN][0]*flter_k[NN][1]+rhok[NN][1]*flter_k[NN][0];

//					rhok[k+nz*(j+ny*i)][0]=rhok[k+nz*(j+ny*i)][0]*flter_k[k+nz*(j+ny*i)][0];
//					rhok[k+nz*(j+ny*i)][1]=rhok[k+nz*(j+ny*i)][0]*flter_k[k+nz*(j+ny*i)][0];

				}
			}
		}
	}
	fftwf_execute(backward);

	r=&phi_out[0][0][0];

	for(i=0;i<nx*ny*nz;i++){
		r[i]=phi[i][0];
	}
	fftwf_free(rho);
	fftwf_free(rhok);
	fftwf_free(phi);
	return 0;
}
